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I. INTRODUCTION 

Interesting processes in physics are frequently associ- 
ated with nonequilibrium situations. For instance, for 
a device to do some work an energy current must flow. 
The state of the system is therefore a nonequilibrium 
one. It is clearly desirable to understand nonequilib- 
rium physics, unfortunately, as opposed to equilibrium 
physics, no general theory of nonequilibrium processes 
exists, for instance, for their stationary distribution. For 
equilibrium systems on the other hand the invariant equi- 
librium distribution is well known. The main difference 
between equilibrium and nonequilibrium systems is that 
in an equilibrium a principle of detailed balance holds, 
stating that the net flow of probability between any two 
states x\ and X2 is zero. That is, in equilibrium the 
probability to go from x\ to X2 is equal to the prob- 
ability to go from xi to x\. This simple rule greatly 
simplifies the analysis. In a stationary nonequilibrium 
situation on the other hand, only the total probability 
flow out of a state x\ has to be zero and not along each 
connection individually. In view of the lack of general 
theory it would be desirable to find some exact solutions 
for nonequilibrium situations, from which we could per- 
haps draw some generic rules. Particularly simple are 
nonequilibrium states that do not change in time, also 
called nonequilibrium stationary states (NESS). 

In classical physics there are a number of exactly 
solvable nonequilibrium models, most notably various 
stochastic lattice gasses described by exclusion pro- 
cesses [IHIj]- The picture is quite different in quantum 
physics. There are hardly any analytically solvable mod- 
els known, see though for instance an example of a single 
spin coupled to a bath in a star configuration An 
exception of solvable models are those that are quadratic 
in the fermionic variables. One possibility for studying 
NESS is to take an infinite system in which an infinite 



sub-part serves as a bath. Such was the case in the dou- 
bly infinite XY chain studied in Q. The other approach 
is to write an effective master equation that describes the 
evolution of only the central system without baths. The 
simplest master equation is the Lindblad equation, which 
can be diagonalized if the superoperator (a generator of a 
markovian flow) is a quadratic function of fermi onic op- 
erators 0] , giving the exact solution for the NESS @, Q ■ 

In the present work we are going to find a solution for 
quantum NESS in a system whose Lindblad superopera- 
tor C (EqJ2|) is not a quadratic function of fermions but 
is, nevertheless, simple enough to enable an explicit so- 
lution. Some parts of the solution have been presented 
in recent Refs. [t| [l(|, while the same model has been 
numerically studied in (lTj . The model we are going to 
study is the XX chain with a local dephasing at each 
site and coupled to nonequilibrium baths. A nonequi- 
librium XX model without dephasing was solved in [12| |. 
for a compact solution see also Besides providing 

an explicit solution we also show that it is fundamen- 
tally different from those in quadratic systems. Namely, 
the Wick theorem does not hold and the NESS is not 
Gaussian. Another interesting aspect is that in a certain 
thermodynamic or weak driving limit the solution can be 
written as a matrix product operator (MPO) with matri 
ces of small fixed dimension 4. This extends the MPO so- 
lution obtained for the NESS in a model without dephas- 
ing ficl ] . The applicability of a matrix product ansatz has 
important implications for the properties of the system 
as well as for numerical methods that can be used to 
calculate NESS. It has been known for some time that 
the matrix product ansatz can describe ground states 
of certain 1-dimensional systems [TH, NESS of classical 
exclusion models 0, [l4| , and is a rather useful concept 
in quantum information and numerical methods used to 
simulate quantum systems [l5[. A MPO ansatz of small 
fixed dimension can describe the time evolution of certain 
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operators in intcgrable systems [161 - 4181 ] . The MPO solu- 
tion presented here extends the applicability of a MPO 
ansatz to NESSs in non-quadratic systems. 

The final section deals with phase transitions in 
nonequilibrium stationary states. It is shown that the 
model exhibits a nonequilibrium phase transition at zero 
dephasing, going from the NESS that is Gaussian and 
displays a ballistic transport, to the NESS that is non- 
Gaussian, shows diffusive spin transport, and exhibits 
long-range correlations for nonzero dephasing. We show 
that the characteristic feature of the phase transition 
point is a faster closing of the gap of the superoperator 
with the system size than away from the phase transition. 
This seems to be a general property of nonequilibrium 
phase transitions, see recent studies in Refs. fA flji]. 



The dephasing part £ dc P h — ^™ =1 £^ cph is a sum of 

C^ eph , each of which acts only on the j-th site and is 
described by a single Lindblad operator, 



- deph 



(5) 



Sometimes it is useful to write a matrix representation 
of the dissipative superoperator £^ eph . If we use a basis 
of Pauli matrices and we order them as {cr*, crj, ct 2 , lj}, 
we have a matrix representation 



idcph 
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II. THE MODEL 

The Hamiltonian of the XX spin chain in a homoge- 
neous magnetic field is given by 



(i) 



i=i 



with standard Pauli matrices. An exact evolution of a 
central system, like the XX chain in our case, which is 
coupled to environment is in general complicated. How- 
ever, if the environment correlation time is sufficiently 
small, i.e., so that it has no "memory", and the evolu- 
tion has a property of the semigroup, then the reduced 
evolution of the central system can be described by the 
Lindblad equation [20j . 



-p = i[p,H] + C dis (p) = C(p). 



(2) 



The dissipative linear operator £ dls can be written in 
terms of Lindblad operators Lfc , 



C dis ( P )=J2(l L kP,Ll} + [L k ,pL{] 



(3) 



The Lindblad master equation can describe the most 
general completely positive trace preserving map that is 
a dynamical semigroup, i.e., a map that is a semigroup 
for a continuous time parameter. While the complete 
positivity must arguably be satisfied by any evolution, 
the semigroup property can be violated when the envi- 
ronment has a "memory" that causes a back-action on 
the central system. 

The model we study here has two different dissipa- 
tive parts. One describee the action of two baths, in- 
ducing a nonequilibrium situation if they are different, 
while the other describes a local dephasing at each site, 
for instance, being due to the coupling of each site to un- 
observablc degrees of freedom. The Lindblad dissipator 
is therefore a sum of two terms, 



-•dis 



bath 



£dcph 



(4) 



Dephasing with strength 7 causes an exponential decay of 
the off-diagonal elements of a density matrix, if we write 
it in the diagonal basis of cr z . Using Jordan- Wigner trans- 
formation, such basis corresponds to the number basis of 
spinless fcrmions. Note that it is precisely the dephas- 
ing term which makes the superoperator non-quadratic, 
in fact quartic, in fermionic operators. Namely, the su- 
peroperator £ de P h is quadratic in the Lindblad operator 
oc cr z which is itself quadratic in fermionic operators. 

For the dissipative bath part £ bath we shall take the 
simplest possible Lindblad operators that are still able 
to describe a nonequilibrium situation. First, they are 
going to act locally only on the first and the last spin, 
and second, there will be only two operators at each end. 
Writing the dissipative part C 



bath 



/•bath 



'•bath 



as 



a sum of a part acting only at the left end (site index 
j = 1 and label "L" ) and a part acting only at the right 
end (site index j = n and label "R" ) , we have (p) = 

E*=i, 2 (l^'V ifc' Rt ] + [£fc' R , ^ ,Rt ]) • The two Lind- 
blad operators at the left end are 



(7) 

while on the right end we have 



Lf = ^T R (l + p + fi)a+, L^= y /r R (l- f i-p)a-, 

(8) 

(Tj = (a? ± i<rJ)/2. The matrix representation of both 
superoperators is 



-•bath 



£bath 
R 




r R 



/ -2 
0-20 
00-4 

V 






-i(-H-p) 




(9) 



(10) 



Such simple local Lindbald operators involving a + and 
a~ are often used when studying transport in spin chains, 
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sec for instance [2lM25| . Two parameters IY and Tr play 
the role of a coupling strength, while p and p determine 
the magnetization that the bath tries to impose on the 
chain ends. To see that this is indeed the case one can 
look for a stationary state of the bath £^ ath dissipator 
only. That is, we look for a single-spin state p that satis- 
fies £k ath (/5) = 0. One easily finds that p ~ li — (p,— p)a\ 
and therefore tr (paf) = —p + p. The chosen bath is 
therefore such that it tries to induce a local magnetiza- 
tion of size ~p + p at the left end and p + p at the right 
end (these two parameters can be thought of as magneti- 
zations of an infinite bath to which the central system is 
coupled) . Of course, the nonequilibrium stationary state 
of the whole master equation ([2]), which in addition in- 
cludes a unitary part and a dephasing, will have a slightly 
different magnetization at the ends. One can also invert 
the problem and ask, is it possible to choose Lindblad 
operators that will target an arbitrary stationary state p 
of £ bath j even on many qubits? The answer is positive 
with the explicit procedure given in [24 1 . 

The goal of this paper is to analytically find the NESS 
of the master equation $Z§. This state, simply denoted 
by p in the following, is an eigenstate of the Lindblad 
superopcrator with an eigenvalue of 0, C(p) = 0. 



III. THE SOLUTION 

As shown in Ref Q, when p = the solution, i.e. the 
NESS, can be sought in the form of a series in powers of 
the driving p. While a pcrturbativc expansion is always 
possible, our solution is different. It is nonperturbative 
due to the special algebraic structure. We will discuss 
this point in more detail later. The ansatz for the NESS 
is 

p = ± (l + pRU + fn?) + --- prRir) + ...). (11) 

As we shall see, the term p r R^ is of the order p r and 
there is a closed set of equations that give p r Ry). The 
solution can therefore be obtained term by term, without 
having to to deal with the whole set of exponentially 
many equations. We are first going to discuss the case 
with the zero bath magnetization offset, p = 0. The non- 
zero case will then be obtained as a simple modification 
of the solution for p = 0. 

It was shown in Ref. @ that the solution p is a sum 
of terms, where each is a product of operators a z and 
jk = 2(a k a k ' +1 — <J y k a^ + i) at different sites, jk is a spin 
current operator. This simple structure is a consequence 
of the fact, that all other operators, that could result in 
jk or a z when operated on by C, arc zero. The NESS 
can therefore be sought within the algebra of a z and jk- 
Before actually going to the solution itself, let us discuss 
the role of a homogeneous magnetic field of strength B. 

The action of a single B a z term in the Hamiltonian is 

simple. The matrix of the superoperator C rj has only 



two nonzero elements. They arc C^\af) = 2B a y , and 
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-2B aj, while all other are zero, Cj (o^' 1 ) 



0. Because the NESS is a sum of products of only a z and 
currents jk, and because for homogeneous field we have 

J2jCf\jk) = as well as £, rfWfc+i) = 0, the 
NESS for zero field, B = 0, is also an exact solution for 
arbitrary nonzero field B ^ 0. In short, the homogeneous 
magnetic field has no influence on the NESS. 



A. Hierarchy of connected correlations 

Let us briefly argue why the NESS can be calculated 
term by term in the expansion over p, Eq. pip , and why 
the set of equations for each term is closed. Without loss 
of generality we assume that p = 0. A key is a simple 
algebra generated by various superoperators in the mas- 
ter equation. The dephasing term acts as C deph (jk) = 
-4 7 j fc and £ dcph (cr z ) = 0, as well as C dcph {l 3 ) = 0. 

The bath at the left end acts as C^ th {ji) = -2ji, 
££ ath (crf) = -4crf, ££ ath (li) = -Ap<j\. Similar ex- 
pressions hold for the bath at the right end. The uni- 
tary part due to H acts as £j? fc+1 (jfc) = 8(er£ +1 — of,), 
£ H K) = (jk - jk-i) and £ H (l fc ) = 0. Note that 
t here are no overlapping products of operators in p. 
The condition that there are no products of operators 
at the same site, i.e., either o z a z or jkjk, ensures that 
the explicit normalization of p by 1/2™ is not broken 
(jlljl . Other terms not present due to hcrmiticity are 
jkcrl + aljk = and jk<r%+i + <J z k+1 j k = 0; note however 



that jkjk+i + jk+ijk 



are allowed. From the action of superoperators we note 
two things: (i) if r is a sum of operators, each of which 
is a product of non-overlapping <r| and jk, then £(r) is 
again a sum of non-overlapping a z and jk', and (ii) the 
number of operators <r| and jk in each term is preserved 
by the action of £ H and £ d °P h . The only superoperator 
that does not conserve the number of operators is that of 
the bath, which can create a z out of an identity. These 
two observations are crucial. Let us now write our solu- 
tion p as a series in p, Eq. where the r-th order term 
is a sum of all possible non-overlapping products, where 
each is a product of exactly r operators a z or jk (here 
we always mean the number of non-identity operators). 
We put all (unknown) coefficients appearing in the term 
pTRir) - m a ge £ £(r) _ -p or ms i ance; the set has only 

one coefficient, namely a known normalization 1/2™ in 
front of 1, the set consists of n unknown coefficients 
in front of a z and n — 1 unknown coefficients in front of 
jk- The NESS must satisfy the equation C(p) = 0, there- 
fore, the coefficient in fron t of each operator appearing 
in C(p) must be zero. If we look at the coefficient in 
front of an operator that is a product of r non-identity 
operators, then this coefficient is a linear function of co- 
efficients in the set S^~' and coefficients in the set S f ' r_1 ' 
(these come from the action of the bath £ bath ). This 
structure enables us to solve for p r R^ r \ i.e., coefficients 
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in S( r >, iteratively, starting with known S^ '. Note that 
the only coefficient in is equal to 1/2" ~ fjP. Using 
this known S 1 - " 1 we can now write a set of linear equa- 
tions for coefficients in where the coefficient from 
will act as a source term, i.e., an inhomogeneous 
part of equations. Because when £ bath makes a term of 
type 5'( r + 1 ) from S( r \ it always multiplies it by //, the 
source term will scale as /i 1 . Therefore, all coefficients 
in S^ 1 * 1 are proportional to fx 1 . Iteratively repeating the 
procedure we see that, (i) coefficients in scale as fx r , 
that is the terms in /j, T R^ indeed scale as [i r , and (ii) 
at each order we have to solve a closed set of equations 
for coefficients in S^- r \ These determine all r-point cor- 
relations in the NESS. Helping ourselves for the moment 
with the solution for the first three orders given in the 
following, we can say even more. If one writes equa- 
tions for r-point connected correlations, instead of for 
non-connected ones, we can also predict how the source 
term at the r-th order scales with the number of spins n. 
Using equations from the previous order S^^ 1 ^ one can 
see that the source term scales as /i r /rt r . In the following 
we shall see that the spin current, the coefficient of which 
is denoted by b, scales in the same way. Therefore, the 
source terms scale as the r-th power of the current, ~ b r . 
A consequence of this is that the largest co nnected term 
in /i r i?( r ' scales as ~ b r n and comes from the connected 
correlation of r operators u|. The largest non-connected 
term in iJ, r R^ on the other hand scales as ~ b r n r and 
comes from the product of r operators see also the 
explicitly ansatz below. 

We are now going to find an explicit form of the first 
three orders. For special values of parameters, the solu- 
tion for the first two orders has been presented in Q. 



B. No magnetization offset, fi — 

Let us first find the NESS for a bath with a zero offset 
of magnetization, fi = 0. It is instructive to first find 
the equilibrium stationary state in the absence of driv- 
ing, when /i = 0. In this case both baths ([7]) act with a + 
and <j~ with equal probability, inducing no net magneti- 
zation. In fact, the equilibrium state is very simple and 
equal to the totally mixed state, 



1 



Poq 



(12) 



This equilibrium state can therefore be thought of as an 
infinite temperature state. 



1. First two orders 



The Ansatz for the first two orders is the following, 

n ^ n— 1 

3=1 k=l 



The spin current operator is j k = 2(a\a\ +l — o\a\ +l ). 
The 2nd order ansatz is 

M 2^(2) _ M_ _|_ BA) + ^2 C + ^2 D + ^a_pi ( 14) 

n n 

H 2 C = Y^ {C 3 , k +a 3 a k )cj]al, 

3=1 k=j+l 



n-2 



n-1 



n—l—j 



f D = E i E <#« - E I > 

3=1 \l=j+l 1=1 

^ F =gYl 3k3l- 

k,l=l 



(15) 



In the solution a specific factor will appear in all expres- 
sions for connected correlation functions. We denote it 
by the letter t, 

. _ (£l + r R )(i + r L r R .) + 2(n - 2) 7 r L r R 

2 7 r L r R ■ [W) 

For large n it scales as t ~ n, in fact, for Tl = Tr = 1 it 
is equal to n — 2 + 2/j. It therefore plays the role of an 
effective system size. 

The solution for the 1st and the 2nd order is found in 
exactly the same way as in Ref. Q for specific parame- 
ters. One writes a closed set of equations and solves it. 
We do not repeat the details of the derivation here, but 
just present the solution for general values of parameters. 
The 1st order terms are 



b = —fx 



and 



2r L r R 



(r L + r R ) (i + r L r R ) + 2(n - 1)7 (t + 1)7 

(17) 



01 



0,2 = ax - b(T L + 2 7 ) 
«3 = «2 - 2 7 6 

Cln-1 = a n-2 - 276 

o« = a n -i - &( r R + 27) = A f + 



Ft 



(18) 



Alternatively, we can express local magnetizations Oj as 
a 3 =-n-b kf ] =n + b kf ] (19) 

fcf ) = {^ i ^ + 27,..., i ^ + 2(n-l)7 + r R }, 



r L ' r L 



Fl 



JR) -{^ + 2(n-l)7 + F L ,...,i±il + 27,^}. 

1 R 1 R 1 R 



(20) 

Away from the boundary we have fey+i — k^ = 27, while 
, (R) _ j,(R) _ 



27. The expectation value of the magneti- 
(13) zation at site j is just (a 7 -) = a 3 and of the spin current 
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{jk) = 26. Several interesting observations can be made. 
Our one-spin Lindblad bath is constructed in such a way 
that it targets magnetization =p/i on the 1st and the last 
site. Because there are other terms besides the bath in 
the Lindblad equation the actual magnetization in the 
NESS at the boundaries is slightly different. Namely, we 
have oi = — b/T^ and a„ = Li+b/T^. Expectedly, one 
can see that if the coupling strength IY (or Tr) is very 
small the difference from the targeted magnetization is 
large. Second, the magnctizaton difference between two 
neighboring sites is 2<yb in the bulk, while it is b(2-f + 
at the left end and b(2j + Tr) at the right end. The 
magnetization profile is theref ore linear, apart from two 
sites at both ends, where a contact resistance due to IY 
and Tr is felt. One could in fact generalize our solution 
to allow for an inhomogeneous dephasing jj at each site. 
The difference in the first order solution (|18|) would be 
having terms "fj + 7j+i instead of 27. Using 71 = j n = 
and having IY = Tr = 7 one could therefore achieve 
a perfectly linear profile, without any contact resistance 
jumps at the boundaries. However, higher order terms in 
/i are more complicated in this case (although of the same 
form) and we do not discuss inhomogeneous dephasing in 
the present work. 

The expectation value of the current 2b is maximal for 
intermediate couplings IY.r. Keeping the system size 
and the dephasing strength 7 fixed, we can see that 
b has a maximum at IY = Tr = f. It is equal to 
\b m ax\ = + (n — 1)7). If the couplings are smaller 
the current is smaller because the NESS is only weakly 
nonequilibrium due to the weak coupling. On the other 
hand, if Ts are large the current is again smaller be- 
cause the hamiltonian part H, which transports mag- 
netization, is less important compared to baths. The 
dependence of the maximal current on 7 is trivial, the 
smaller the dephasing the larger is the current. Another 
observation is that the magnetization difference between 
the two ends is a n — a\ = — fo(IY + Tr + 27(71 — I)). 
The transport coefficient k, defined via j = — K a "~ ai , is 
K = 2n/(IY + Tr + 2{n — 1)7). As long as the dephasing 
is nonzero it asymptotically behaves as k ~ I/7 and is 
independent of the system size n. If 717 is on the other 
hand smaller than IYjIY, for instance if the dephasing 
is zero, then in the limit of weak coupling, IY,rR — > 0, 
k diverges. 

The expressions for the 2nd order terms are rather sim- 
ple, 

/ = 6 2 (1 + 1/t) (21) 

d ^ { ; j > 1 (22) 

T \ fc™ ; j <i (22) 

Cij = -j (fc^fcf 5 + (1 + t) S i+hj ) . (23) 

The connected correlation function of magnetization, 
(cr-cr-)c = (of<rf) - (<J-)(<Jj) is for i = j equal to I - af, 



while it is equal to (of o|) c = Cjj for nondiagonal i ^ j. 
For large system size n the connected correlation func- 
tion Cij achieves its maximal value at i f» j w n/2, when 

a product of k^kj is the largest. We therefore have 
max(Ci.j) ~ ( r )nb) 2 /t. The maximum is again achieved 
at Tl = Tr = 1. The dependence of the maximal z — z 
connected correlation on 7 is though the opposite of that 
of the current. Here the correlations monotonously in- 
crease with increasing 7. 

The 3rd order terms are obtained in an analogous way. 
The calculation is a little tedious, nevertheless, we man- 
aged to obtained closed expressions. Detailed results can 
be found in the Appendix „ At this point we just list 
the dominant term for large n. It is a three-point con- 
nected correlation of o z . Away from boundaries and 
for non-neighboring indices, where the Kronecker-delta 
terms are zero, the general expression (|A.4[) simplifies. 
For instance, for 7 = Tl = Tr = I and i < j < k it is 
just 

(o*o* j0 %) c = -^-rx(l-2y)(l-z), (24) 
J n[n — 1) 

where we have introduced rescaled position variables x = 
—^r , y = -^rr , z = — tt ■ The two-point function is for the 
same parameter values and i < j 

(*l*°) c = -^x(l-y). (25) 
J n 

Scaling of the first two connected correlations on \i and 
n therefore follows a general rule: the r-point connected 
correlation function is ~ [i T /n r ~ 1 in accordance with the 
general discussion in the subsection 1111 Al The form of 
these two correlation functions is the same as in some 
classical exclusion processes Q, for a recent work on an 
analogous quantum master equation exhibiting some fea- 
tures similar to our model see [2?} ■ Whether this is only a 
consequence of the same hydrodynamic limit or whether 
there perhaps exists an exact mapping from the quantum 
model to a classical one is unknown at present. 

C. Nonzero offset, jl 7^ 

Let us now go to the case where there is an offset of 
magnetization in the baths. For nonzero fx the equilib- 
rium state, when \x = 0, is again simple, although not 
trivial ~ 1 as in the case of zero jl. It is a product state 
with nonzero magnetization 

1 " 

3=1 

It is therefore a state with all connected correlations be- 
ing zero (trivially Gaussian state) apart from nonzero 
magnetization. This equilibrium state can be equated to 
the grandcanonical state p gc ~ exp (—(3(H — x^ z )); with 
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£ z = ■ (Xj , having an infinite temperature ft = and a 
finite chemical potential times the inverse temperature, 
0X = \ l n ((1 + — £*))■ One should note though 

that integrable systems, such as our XX chain, do not 
thermalize for generic local Lindblad baths [26| . 

Out of equilibrium, when [i ^ 0, the solution is actually 
very similar to the one with zero p,. Namely, the only 
thing that changes is the expression for magnetization, 
a,j , while all other n-point connected correlations are the 
same as in the case of no offset, p, = 0. For nonzero offset 
the solution is therefore 

dj =p-n-bkf J \ (27) 

with the same k^ (|2"U| as before. All other terms 
(|17I23IA.4IA.5IA.6P are the same. 

IV. THE NESS IS NON-GAUSSIAN 

In this section we are going to show that the NESS ob- 
tained above is not Gaussian, i.e., that the Wick theorem 
does not hold. Because three-point connected correlation 
functions are nonzero, see e.g. eq. (|A.4[) . it is obvious that 
the Wick theorem does not apply in spin variables. How- 
ever, it is not clear that it does not apply in the spinlcss 
fermion picture. 

A system of spin-1/2 particles can be mapped to 
spinlcss fcrmions using the Jordan- Wigner transforma- 
tion. Denoting by Cj and canonical fermionic annihila- 
tion and creation operators, satisfying anticommutators 
{cj,Cfc} = 0, {c],c\} = 0, {cj,c\} = Sj.kl, the transfor- 
mation is given by the mapping 

Cj = -(^•••<7)_>/ 

4 = -K---a]_ 1 )aJ, (28) 

or its inverse 

*i = -(*;• "Oj-ijfe+ct) 

^ = -i(af-^_ 1 )(c J --c}) 

erf = Cjc] - c]cj = 1 - 2n h (29) 

where we denote by rij = CjCj a number (density) oper- 

(r) 

ator at site j. Denoting by Z - ' = a % y ■ a^ +r _ 1 a string 
of r consecutive er z , and introducing an energy-density- 
like operator H<f +1) = a* aj +r + a] Z^ aj +r 

and a current-like operator B^ +1 ^ = aj Z^ + ^a^ +r — 
ffji?^ cr^,, we have 

H ( f^=2{c]c j+r -c j c) +r ) 
i?j r+1) = 2i{c]c J+r + c 3 c] +r ) 

HJ^=2{c]c) +r -c j c j+r ) 

BJ {r+l) = 2\{c]c] +r + c jCj+r ), (30) 



where H~ (r+1) = a*Z^o* +r - a^Z^a y j+r and 

BJ [r+1) = a^zf+^a]^ + ojzj+i 15 ^. From our so- 
lution for the NESS (see also comments in ref. (@)) 

we can see that expectations of all Hj T \Bj T+1 \Hj 

— (r) (2) 
and Bj , apart from B - , are zero. The only nonzero 

two-point fermionic expectations are therefore (CjCj) = 

(1 — aj)/2 and 2i(c]c J+ i + Cjcj +1 ) = b. An important 
point is that only on-site or nearest-neighbor two-point 
fermionic correlations are nonzero. If the NESS would 
be gaussian and the Wick theorem would hold all con- 
nected correlations beyond nearest neighbor would have 
to be zero. Because this is not the case the NESS is 
clearly not gaussian. For instance, rewriting the con- 
nected z — z correlation as Cij = (c|cr|) — (er z )(er z } = 

4((c]ci CjCj) — (c-Cj)(cjcj)), and using the Wick theorem, 

we would have Cij = — (cjcj)(cjct). The last expression 
is nonzero only if j = i+ 1, whereas on the other hand we 
have long-range correlations with all Cjj being nonzero 
if 7 7^ 0. The Wick theorem therefore does not hold if 
7^0. The NESS is non-Gaussian and presents an in- 
teresting new solvable model. Whether it is nevertheless 
equivalent to some existing solvable model is at present 
unknown. 



A. Matrix product operator ansatz 

Even-though the NESS is non-Gaussian, it is still, in 
a sense, weakly correlated, meaning that it can be rep- 
resented in terms of a matrix product operator (MPO) 
form 

p = ^ zZ^iA^A^ ■ ■ ■ |i> ■ ■ ■ < n , 

(31) 

with matrices of small size. Note that in an MPO for- 
mulation a density matrix is treated as an element (a 
pure state) of 4" dimensional Hilbcrt space of opera- 
tors. For 7 = the matrices A" of finite size D = 4 
that is independent of the system length suffice [l(| • For 
nonzero dephasing an exact representation with matrices 
of fixed size is not possible. We have numerical indica- 
tions though that the Schmidt rank, i.e., a number of 
nonzero Schmidt coefficients, for a bipartite cut after the 
first m spins is 4m. We have verified this conjecture by 
numerically computing an exact NESS for systems of up- 
to n = 10 spins. Because the Schmidt rank is equal to 
the necessary dimension of matrices in the MPO ansatz, 
the exact representation of the NESS requires an MPO 
of dimension D = 2n. 

Looking at the series solution for NESS (fTTj) we can see 
that the largest connected term in the r-th order fj, r R^ 
(i.e., r-point connected correlation function) scales for 
large n as ~ b r n = l/n r ~ 1 and comes from the con- 
nected correlation of r er z s. Therefore, in the thermo- 
dynamic limit of large n one could neglect all higher or- 
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der connected correlations, keeping in the NESS only the 
leading terms that scale with n as ~ 1 and ~ 1/n. These 
are magnetization <r| , spin current jk and z — z connected 
correlations, as well as their products decaying no faster 
than ~ 1/n that appear in higher order non-connected 
correlations. Because all these terms are already present 
in an exact MPO solution for 7 = having D = 4, 
we can expect that in the limit of a large system one 
might be able to approximately describe the NESS with 
an MPO of dimension smaller than 2n. An MPO with 
dime nsion D = 4 that correctly describes all terms larger 
than or equal to ~ 0(1/ n) in the NESS is a simple ex- 
tension of the solution for 7 = [lOj . Four matrices at 
each site have to have the following form, 



H=0.1 




-4 



(i) 



sin -7) P, A 
2 JI 



(y) 
3 



P = 



where 



/OflO 



1 

Vo 

- —k 







(L) 



R = 



( 1 
10 


Vo 

= -(c°s 2 3 + sm -^3) 

4 0' 





10 



R 



(32) 



and 



4fc 



R) 



Parameters a. 



Matrices A- x) and A (y) 



.4 



(x) 



i+4 — ^ , and can be 
(-P,-P,P,P,-P,...), and 
. .), if one writes them as a 



and b are given in Eqs. (|18ll7[) . 
are periodic with period 4 
concisely written as = 
A&> = (-R,R,R,-R,-R, 
vector. 

Whether the above MPO with D = 4 actually suf- 
fices to describe the NESS in the thermodynamic limit 
depends on the scaling of the next largest Schmidt coeffi- 
cient A5 (in this subsection we denote by Xj, j = 1, . . . the 
Schmidt coefficients for a symmetric bipartition at n/2 
spins; J2j = !)• If -D = 4 is to suffice, A5 must decay 
more rapidly than A4. Because the scaling of Schmidt co- 
efficients is not simply related to the scaling of expansion 
coefficients of the NESS we used numerical simulation to 
study the scaling of Xj. We have calculated the NESS 
using a time dependent density matrix renormalization 
group method [28| [TTJ] (tDMRG) for systems with up-to 
n = 128 spins using an MPO ansatz with a fixed small 
matrix dimension D = 10. In all results that follow we 
have fixed Tl = Tr =7 = 1 and p, = 0. As an indepen- 
dent check of our asymptotic MPO (|3"2"j) solution, we also 
compared the f our largest Schmidt coefficients obtained 
from the tDMRG with those obtained from the D = 4 
MPO (|32|) . Having an explicit representation of matrices 
(f3"2"|) it is easy to calculate the Schmidt coefficients. 

The results for Xi t ,„fi(n) are shown in Fig. [TJ Sev- 
eral interesting things can be observed. First, the largest 
four Schmidt coefficients obtained by the tDMRG agree 
with the analytical calculation from Eq. (|3"2"|) . Interest- 
ingly, around n rj 60 for a = 0.1 and around n w 30 for 
fx = 0.2 an avoided crossing occurs between A2,3 and A4. 
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FIG. 1. Scaling of the 6 largest Schmidt coefficients Xi of the 
NESS with the system size. Top figure is for /x = 0.1, bottom 
for (i — 0.2. Symbols are tDMRG results for Ai, A2, A3, A4, A5 
and A6 (stars, open circles, squares, crosses, triangles and 
full circles, respectively); full lines are analytic Schmidt co- 
efficients from MPO with D = 4 (|32J while the dotted line 
indicates asymptotic ~ 1/n scaling of A5. We use a symmetric 
bipartite cut after n/2 spins, IY = Tr = 7 = 1, p, = 0. 



A consequence of this is that for large n we have scaling 
A4 ~ 1/n, while we have A4 rj n° for small n. From 
the data for different /z we can infer that the crossing 
happens when an = v c , with v c w 6. Similar avoided 
crossing at the same value of n occurs also between A5 
and Ag . For an < v c A5 decays faster than 1 /n while it 
decays as A5 ~ 1/n for (in > v c . We therefore see that 
if we want A5 to decay with n faster than A4, in other 
words, for D — 4 MPO (|3"2"|) to really give the leading 
term solution, one must have fin < v c when going to the 
thermodynamic limit n — > 00. 

Another possibility for the D — 4 solution to give the 
leading order would be a limit of small driving, fi — > 0. 
If A5 decays with fi faster than A4 this would still be 
enough. We therefore also looked at the scaling of Aj with 
fi. Looking at Fig.[5J and noting that the transition point 
of the avoided crossing fin = v c happens at fi m 0.05 for 
n = 128 and fi « 0.1 for n = 64, we can observe the 
following scaling for the largest coefficients: (i) A2 ~ fi 

p; 



for fin < v c , while A2 ~ fi 2 for fin > v c ; (ii) A3 
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FIG. 2. Scaling of the Schmidt coefficients A2.3,4.5 with /j,, 
obtained by the tDMRG (circles) for n = 8, 16,' 32, 64, 128. 
Dashed lines suggest the asymptotic scaling with /i. We use 
a symmetric bipartite cut after n/2 spins, IY = Tr = 7 = 1, 
p, = 0. 



(iii) A4 ~ /z 2 for fin < is c , while A4 ~ fi for fin> v c ; (iv) 
A5 ~ /i 2 for /in < iv c , while A5 ~ /1 3 for un > v c . We 
can see that for fi n > i/ c A5 decays faster with u than 
A 4 . The MPO with D = 4 J32]) therefore also gives the 
leading order solution in the limit of weak driving, u — > 0, 
provided we have fin > v c . Presumably there are further 
avoided crossings in the spectrum of Schmidt coefficients, 
besides the one at fi n = v c . fin should therefore be 
smaller than the value at these higher crossings, however, 
these points need further investigation. 

To summarize, an MPO ansatz of size D = 4 (pi!?]) gives 
the leading order solution in two limits, either n — > 00 
while keeping fi n < v c , or fi — > while keeping fin > v c . 
The two conditions can in fact be put under the same 
hood. Observing that to fulfill nfi < v c in the ther- 
modynamic limit one must necessarily have fl — Y 0, and 
similarly, to have n \i > v c in the weak-driving limit one 
must necessarily have n — > 00, one can reformulate both 
in a single statement, saying that the MPO of dimension 
D — 4 gives the leading order solution in the limit 
n — » 00 having at the same time \i — > 0. 

Another interesting point is that for fin > v c the sec- 
ond largest Schmidt coefficient is independent of n and 
scales as A2 ~ /i 2 /n°. This means that despite the fact 
that our system is at an infinite effective temperature, 
there is a nonzero bipartite entanglement present, even 
in the thermodynamic limit of large n (at a fixed fi). 
This entanglement at an infinite temperature is of purely 
noncquilibrium origin. 



V. NONEQUILIBRIUM PHASE TRANSITION 

From the exact solution we can see that the NESS 
undergoes a transition from a state without long-range 
correlations for 7 = 0, to the one with long-range corre- 
lations for 7 ^ 0. A point in a parameter space where 
a system undergoes a sudden change in some expecta- 
tion values is usually called a phase transition point. We 
use the same nomenclature here and call this transition 
a nonequilibrium phase transition, because the nature 
of the correlations changes. Transport properties also 
change suddenly at 7 = 0, going from a ballistic (su- 
perconducting) state to a diffusive for nonzero dephas- 
ing. In the phase with long-range correlations two-point 
z — z correlations scale as ~ /i 2 /n Q and are therefore 
of purely noncquilibrium origin. They also go to zero in 
the thermodynamic limit making this transition different 
than equilibrium phase transitions. The correlation func- 
tion has a plateau because in the thermodynamic limit 
the decay of Cij with the distance between indices — 
gets increasingly slower (|25p . Recently, similar quantum 
noncquilibrium phase transitions have been discovered in 
the XY model 0, Q , in which the long-range correlations 
scale as ~ fj 2 /n, in the XXZ model with dephasing [ll| 
(the same scaling of the correlation plateau as here), as 
well as in the XXZ model without dephasing where 
though the correlation plateau scales as ~ /x 2 /n° (the 
plateau is independent of n!) at an infinite temperature 
and as ~ fj 2 jn a at a finite temperature. It has been con- 
jectured that long-range correlations are a generic feature 
of quantum noncquilibrium steady states [19( . 

It would be nice to understand this phase transition 
more in detail. At equilibrium, phase transitions are con- 
nected with the nonanaliticity of the free energy, or cquiv- 
alcntly, zeros of the partition function. Difficulty with the 
noncquilibrium situation is that there is no general the- 
ory, in particual, free-energy-like quantity whose analytic 
property could be studied is not known. Nevertheless, 
because observable properties of the system change sud- 
denly at our nonequilibrium phase transition point there 
must be some underlying nonanalytic property. We can 
mention that the Lee- Yang theory of equilibrium phase 
transitions has been used with some success on certain 
classical nonequilibrium systems [H, H3, HH whose sta- 
tionary state can be represented in a matrix product 
form. For some classical exclusion processes in the ther- 
modynamic limit even a free-ener gy l ike functional can 
be calculated 3l[. Until recently no quantum sta- 
tionar y state solvable by matrix product ansatz has been 
known and the theory of quantum noncquilibrium phase 
transitions is lagging behind the classical. For classical 
nonequilibrium systems a product of all nonzero eigen- 
values of a superoperator can play a role of a "partition" 
function, Z = EL ^o(~\/)' wnose zeros then determine 
the location of nonequilibrium phase transitions. It is 
shown (33 | that the product of nonzero eigenvalues is 
equal to a sum of all expansion coefficients of the NESS. 
There are problems with such an approach though. The 
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principal difficulty is that Z is defined only up-to to an ar- 
bitrary normalization factor and it is a priori not known 
which normalization one should take. 

In the present work we look for non-analytic signa- 
tures in the spectrum of the supcropcrator, trying to see 
if there is any characteristic change at the phase tran- 
sition point. Because the NESS is an eigenstate of C 
with an eigenvalue zero, a change in the properties of 
the NESS should be connected with the closing of the 
gap between the second largest eigenvalue and the one of 
the NESS. If this gap goes to zero there is the possibility 
for a scenario reminiscent of an avoided-crossing between 
two eigenencrgics happening at a zero-temperature quan- 
tum phase transition, where at the crossing the nature 
of the two levels involved is exchanged. Of interest for 
nonequilibrium phase transitions is therefore the second 
largest eigenvalue A2 of the Lindblad superoperator C, 
or in particular the gap A, A = — A2. We note that the 
size of the gap also determines the decay rate, i.e., the 
relaxation rate, to the NESS. Namely, the initial nons 
tationary state will, for large times, relax to the NESS 
as exp(— At). Because we in general expect the relax- 
ation time, in the case of local coupling to baths at chain 
ends, to grow with the system size at least as oc n - 
larger system simply needs more time to relax -, the 
gap is expected to decrease with n, even if we are not 
at the nonequilibrium phase transition point. The sig- 
nature of nonequilibrium phase transition can therefore 
not be simply the closing of the gap A. What is gen- 
erally observed though is that the closing of the gap A 
with n is at the phase transition point faster than away 
from the phase transition point. Such behavior has been 
observed, for instance, in an open XY chain Q , where at 
a nonequilibrium phase transition point the gap scales as 
A ~ 1/n 5 , while it scales only as A ~ 1/n 3 away from 
the transition point. We therefore studied the scaling of 
the gap in our model. Th e results can be seen in Fig. [3] 
We numerically exactly calculated the gap A for systems 
of up-to n = 11 spins. For larger systems a relaxation 
time r of some observable can be used as an estimate of 
the gap. Using the tDMRG we have solved the Lindblad 
equation in time, thereby obtaining time dependent ex- 
pectations of observables. We fitted an exponential func- 
tion to the relaxation of the magnetization at site n/2, 
z n/2{t) — z n/2(°°) ~ ex P (~ ^/ r )i an d determined r for 
systems with up-to n = 128 spins. As the second largest 
eigenvalue A2 is nondegenerate for 7 = 1 the relaxation 
time can serve to estimate the gap through A ~ 1/r. 
From the figure we can see that at the phase transi- 
tion point the gap scales as A ~ 1/n 3 . This of course 
agrees with a previous analytic result for the isotropic 
XY chain Q . Away from the nonequilibrium phase tran- 
sition, at 7 = 1, the gap decays more slowly, as A ~ 1/n 2 . 
We therefor e conjecture that the characteristic feature 
of nonequilibrium phase transitions is that the gap of the 
superoperator decays faster at the nonequilibrium phase 
transition point than in its neighborhood. At a nonequi- 
librium phase transition point the relaxation slows down, 
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FIG. 3. Scaling of the gap of the Lindblad superoperator 
with system size for the model without long-range order at 
7 = and for the system with long-range order in the case 
of 7 = 1. At the nonequilibrium phase transition point at 
7 = the gap decays faster, as A ~ 1/n 3 , than away from 
the transition point (7 = 1) where A ~ 1/n 2 . We also show 
the scaling of the inverse relaxation time of the magnetization 
at site n/2 (squares). fL = Tr = 1, p, = and u = 0.1. 



similary as at an equilibrium phase transition point. How 
much faster the gap decays depends on a particular sys- 
tem. Here we have 1/n 3 vs. 1/n 2 , for the XY model 
on the other hand one has 1/n 5 vs. 1/n 3 Q. Note that 
this criterion is different from those used in equilibrium 
physics as it makes a reference also to the neighborhood 
of the transition point. Interesting to note is that the 
relaxation time of the spin current is in our model by 
about a factor « 4 shorter than that of z n / 2 - The reason 
must lie in the fact that some of the eigenstates that are 
just below the NESS in the spectrum carry no current. 
Another observation about the spectrum of the Lindblad 
superopera tor for the XX model with dephasing is that 
the eigenvalues are independent of the driving /i, only the 
eigenvectors depend on [i. This probably happens due to 
the hierarchical structure of the master equation. 



VI. SUMMARY 

We have provided exact expressions for all one-point, 
two-point and three-point connected correlations in the 
nonequilibrium stationary state of the XX model with 
dephasing. A hierarchical structure of stationary equa- 
tions is explained, a consequence of which is that equa- 
tions determining all n-point correlations form a closed 
set . The nonequilibrium stationary state is non-Gaussian 
because the Wick theorem does not apply. In the ther- 
modynamic and weak-driving limit the solution can be 
written in terms of a matrix product ansatz with matri- 
ces of fixed dimension 4. At zero dephasing the model 
exhibits a nonequilibrium phase transition from a state 
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with only nearest-neighbor correlations to a state pos- 
sessing long-range correlations. It is conjectured that at 
a quantum noncquilibrium phase transition point the gap 
of the superoperator closes with the system size more 
rapidly than in the vicinity of the transition point. 
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Appendix: 3rd order terms 

We are mainly interested in the connected three- 
point correlations, defined for arbitrary operators A, £?, C 
as (ABC) C = (ABC) - (A) C (B) C (C) C - (AB) C (C) C - 
(AC) C (B) C - (BC) C (A) C . The 3rd order ansatz /i 3 R^ 
is the sum of 5 terms, 



/'•"'/•' : '' = G zzz + G zzj + G zjj + G m + G' m . (A.l) 
The terms are 



n n n 



G zzz = 2J (^*»i» fc + a i a 3°>k + OiCj-jfc + a>jCi t k + a kCij)afafal 

i—l k—j-\-l 

71—1 71 71 - 



G 



k=l i=1 i = i + 1 

i({k,k + l) j({k,k + l} 



; iijfe + 3k ji 



i,k = l 3=1 



E 



:JiJkji 



i,k,l=l 

i^k,ijti,k^i;&t most one overlap 

n— 3 



+3J 



(A.2) 



I 

In the above expressions the summations are such that sites, like jkjk+i- The solutions for unknown coefficients 
no two operators appear at the same site and in addition are 

there are no overlapping terms like jk^k+i- ^ n * nc Cjjj 2 2 

we have at most one overlapping currents on neighboring 

6t(* — 1) ' y *(*-!)' 



Zij,/ — — & 



- {t\kf - + (t + l)(*f ><W - *<%_!)} . (A.4) 



t(f-i) 



fc W (fc (R)_ fc W ) + (t + 1)w 

-fcj R )(fcl R) -fcl L) ) + (t+l)^ +1 

2fc[ L ¥ R) 



k>j 
k < i 

i < k < j — 1 



(A.5) 



*(<-!) 



-2fc 



.(H) 



2A- 



(L) 



L) 



j > fc + 1 

j < i 

i + 1 < j < k 

(A.6) 



Three-point connected correlation function is for nondi- 
agonal indices equal to (afOj<r|) c = Zi,j,k- 
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